LoopDetect - Python¶
Testing in Python 3.12 environment
In [1]:
# Import LoopDetect core and numpy
import numpy as np
import loopdetect.core as ld
import loopdetect.examples as lde
import numdifftools
import pandas as pd
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
d:\Work\StudentJob\Jana Wolf\LoopDetect_2025\env\lib\loopdetect_python_3.12\Lib\site-packages\loopdetect\examples.py:1: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81. import pkg_resources
In [2]:
# An example ODE system with function func_POSm4 is provided in
# loopdetect.examples, it has 4 variables.
# Variable values for the 4 variables (as tuple, cast into a list)
s_star = [(1,1,1,1)]
# Define further arguments of func_POSm4
klin = np.ones((8,))
knonlin = (2.5,3)
# compute loops
res_tab = ld.find_loops_vset(lde.func_POSm4,vset=s_star,klin=klin,
knonlin=knonlin,max_num_loops=10)
# The loop list, a pandas dataframe, is accessed like this:
res_tab['loop_rep'][0]
Out[2]:
| loop | length | sign | |
|---|---|---|---|
| 0 | (0, 1, 2, 3, 0) | 4 | -1 |
| 1 | (1, 2, 3, 1) | 3 | 1 |
| 2 | (0, 0) | 1 | -1 |
| 3 | (1, 1) | 1 | -1 |
| 4 | (2, 2) | 1 | -1 |
| 5 | (3, 3) | 1 | -1 |
In [3]:
type(res_tab) # dict, not a DataFrame
res_tab.keys()
Out[3]:
dict_keys(['jac_rep', 'jac_rep_index', 'loop_rep', 'loop_rep_index'])
In [ ]:
# Don't worry about the following code's errors, it's just for demonstration of some outdated code on our library tutorial.
In [ ]:
klin = [165,0.044,0.27,550,5000,78,4.4,5.1]
knonlin = [0.3,2]
j_matrix = numdifftools.Jacobian(lde.func_POSm4_comp,method="complex")(s_star[0],
klin=klin,knonlin=knonlin).real
signed_jacobian = np.sign(j_matrix)
In [ ]:
# Import the required packages.
import numpy as np
import numdifftools
import loopdetect.examples as lde
# Define kinetic parameters of the model that are arguments of the
# example function func_POSm4_comp provided in LoopDetect
klin = tuple(165,0.044,0.27,550,5000,78,4.4,5.1)
knonlin = tuple(0.3,2)
# Note that we defined s_star in the section above. It could also
# be defined as simple tuple, s_star=(1,2,3,4).
j_matrix = numdifftools.Jacobian(lde.func_POSm4_comp,method="complex")(s_star,
klin=klin,knonlin=knonlin).real
signed_jacobian = np.sign(j_matrix)
In [4]:
# Sort by length of loops found
pd.DataFrame(res_tab['loop_rep'][0]).sort_values(by=['length','loop'])
Out[4]:
| loop | length | sign | |
|---|---|---|---|
| 2 | (0, 0) | 1 | -1 |
| 3 | (1, 1) | 1 | -1 |
| 4 | (2, 2) | 1 | -1 |
| 5 | (3, 3) | 1 | -1 |
| 1 | (1, 2, 3, 1) | 3 | 1 |
| 0 | (0, 1, 2, 3, 0) | 4 | -1 |
In [5]:
print(res_tab['loop_rep'][0])
loop length sign 0 (0, 1, 2, 3, 0) 4 -1 1 (1, 2, 3, 1) 3 1 2 (0, 0) 1 -1 3 (1, 1) 1 -1 4 (2, 2) 1 -1 5 (3, 3) 1 -1
Note:¶
- The
lde.func_POSm4is a network consists of 4 nodes, describe by 4 differential equations - Detaid:
def func_POSm4(x,klin,knonlin):
...
dx = np.zeros(4)
dx[0] = klin[0]-(klin[1]*(1 + x[3]/pow(knonlin[0],knonlin[1])) + klin[2])*x[0]
dx[1] = klin[1]*(1 + x[3]/pow(knonlin[0],knonlin[1]))*x[0] - (klin[3] + klin[4])*x[1]
dx[2] = klin[3]*x[1] - (klin[5] + klin[6])*x[2]
dx[3] = klin[5]*x[2] - klin[7]*x[3]
return(dx)
In [89]:
help(pow)
Help on built-in function pow in module builtins:
pow(base, exp, mod=None)
Equivalent to base**exp with 2 arguments or base**exp % mod with 3 arguments
Some types, such as ints, are able to use a more efficient algorithm when
invoked using the three argument form.
In [5]:
def func_POSm4_fix(x, klin, knonlin):
fb = 1 + (x[3] / knonlin[0])**knonlin[1]
dx = np.zeros(4)
dx[0] = klin[0] - klin[1] * x[0] * fb - klin[2] * x[0]
dx[1] = klin[1] * x[0] * fb - klin[3] * x[1] - klin[4] * x[1]
dx[2] = klin[3] * x[1] - klin[5] * x[2] - klin[6] * x[2]
dx[3] = klin[5] * x[2] - klin[7] * x[3]
return dx
def solve_POSm4(initial_conditions, time_points, klin, knonlin):
def model(x, t):
return func_POSm4_fix(x, klin, knonlin)
sol = odeint(model, initial_conditions, time_points)
return sol
def func_POSm4_comp_fix(x, klin, knonlin):
fb = 1 + (x[3] / knonlin[0])**int(knonlin[1])
dx = np.zeros(4, dtype='complex')
dx[0] = klin[0] - klin[1] * x[0] * fb - klin[2] * x[0]
dx[1] = klin[1] * x[0] * fb - klin[3] * x[1] - klin[4] * x[1]
dx[2] = klin[3] * x[1] - klin[5] * x[2] - klin[6] * x[2]
dx[3] = klin[5] * x[2] - klin[7] * x[3]
return dx
def solve_POSm4_comp(initial_conditions, time_points, klin, knonlin):
def model(x, t):
return func_POSm4_comp_fix(x, klin, knonlin)
sol = odeint(model, initial_conditions, time_points)
return sol
In [6]:
res_tab = ld.find_loops_vset(func_POSm4_fix,vset=s_star,klin=klin,
knonlin=knonlin,max_num_loops=10)
# The loop list, a pandas dataframe, is accessed like this:
res_tab['loop_rep'][0]
Out[6]:
| loop | length | sign | |
|---|---|---|---|
| 0 | (0, 1, 2, 3, 0) | 4 | -1 |
| 1 | (1, 2, 3, 1) | 3 | 1 |
| 2 | (0, 0) | 1 | -1 |
| 3 | (1, 1) | 1 | -1 |
| 4 | (2, 2) | 1 | -1 |
| 5 | (3, 3) | 1 | -1 |
In [7]:
id_to_subject = {
0: "S1",
1: "S2",
2: "S3",
3: "S4"
}
df = pd.DataFrame(res_tab['loop_rep'][0]).sort_values(by=['length','loop'])
df["loop_subject"] = df["loop"].apply(
lambda t: tuple(id_to_subject[i] for i in t)
)
df
Out[7]:
| loop | length | sign | loop_subject | |
|---|---|---|---|---|
| 2 | (0, 0) | 1 | -1 | (S1, S1) |
| 3 | (1, 1) | 1 | -1 | (S2, S2) |
| 4 | (2, 2) | 1 | -1 | (S3, S3) |
| 5 | (3, 3) | 1 | -1 | (S4, S4) |
| 1 | (1, 2, 3, 1) | 3 | 1 | (S2, S3, S4, S2) |
| 0 | (0, 1, 2, 3, 0) | 4 | -1 | (S1, S2, S3, S4, S1) |
In [8]:
# Jacobian calculation
jac_func = numdifftools.Jacobian(lambda x: func_POSm4_fix(x, klin, knonlin))
J = jac_func(s_star[0])
print(J)
signs_J = np.sign(J)
print(signs_J)
[[-2.064 0. 0. -0.192] [ 1.064 -2. 0. 0.192] [ 0. 1. -2. 0. ] [ 0. 0. 1. -1. ]] [[-1. 0. 0. -1.] [ 1. -1. 0. 1.] [ 0. 1. -1. 0.] [ 0. 0. 1. -1.]]
In [9]:
# # Jacobian calculation
# jac_func = numdifftools.Jacobian(lambda x: func_POSm4_fix(x, klin, knonlin), method='complex')
# J = jac_func(s_star[0]).real
# print(J)
# signs_J = np.sign(J)
# print(signs_J)
In [10]:
# Add regulation information for each variable (-1: inhibition, 1:activation). 0 means no regulation and will not be shown in the graph.
# Information is given in Jacobian matrix.
# Try to create a table from the signed Jacobian matrix
# Print to test
for i in range(signs_J.shape[0]):
for j in range(signs_J.shape[1]):
print(f"Regulation from variable {j} to variable {i}: {signs_J[i,j]}")
Regulation from variable 0 to variable 0: -1.0 Regulation from variable 1 to variable 0: 0.0 Regulation from variable 2 to variable 0: 0.0 Regulation from variable 3 to variable 0: -1.0 Regulation from variable 0 to variable 1: 1.0 Regulation from variable 1 to variable 1: -1.0 Regulation from variable 2 to variable 1: 0.0 Regulation from variable 3 to variable 1: 1.0 Regulation from variable 0 to variable 2: 0.0 Regulation from variable 1 to variable 2: 1.0 Regulation from variable 2 to variable 2: -1.0 Regulation from variable 3 to variable 2: 0.0 Regulation from variable 0 to variable 3: 0.0 Regulation from variable 1 to variable 3: 0.0 Regulation from variable 2 to variable 3: 1.0 Regulation from variable 3 to variable 3: -1.0
In [11]:
# Dataframe with columns: 'from', 'to', 'regulation'
# Apply id to subject mapping
regulation_info = pd.DataFrame(signs_J, columns=[id_to_subject[i] for i in range(signs_J.shape[1])], index=[id_to_subject[i] for i in range(signs_J.shape[0])])
regulation_info = regulation_info.stack().reset_index()
regulation_info.columns = ['from', 'to', 'regulation']
# Add regulation information for each variable (-1: inhibition (negative), 1:activation (posotive)). 0 means no regulation and will not be shown in the graph.
regulation_info['regulation_words'] = regulation_info['regulation'].map({-1: 'inhibition (negative)', 1: 'activation (positive)', 0: 'no regulation'})
regulation_info
Out[11]:
| from | to | regulation | regulation_words | |
|---|---|---|---|---|
| 0 | S1 | S1 | -1.0 | inhibition (negative) |
| 1 | S1 | S2 | 0.0 | no regulation |
| 2 | S1 | S3 | 0.0 | no regulation |
| 3 | S1 | S4 | -1.0 | inhibition (negative) |
| 4 | S2 | S1 | 1.0 | activation (positive) |
| 5 | S2 | S2 | -1.0 | inhibition (negative) |
| 6 | S2 | S3 | 0.0 | no regulation |
| 7 | S2 | S4 | 1.0 | activation (positive) |
| 8 | S3 | S1 | 0.0 | no regulation |
| 9 | S3 | S2 | 1.0 | activation (positive) |
| 10 | S3 | S3 | -1.0 | inhibition (negative) |
| 11 | S3 | S4 | 0.0 | no regulation |
| 12 | S4 | S1 | 0.0 | no regulation |
| 13 | S4 | S2 | 0.0 | no regulation |
| 14 | S4 | S3 | 1.0 | activation (positive) |
| 15 | S4 | S4 | -1.0 | inhibition (negative) |
In [3]:
import matplotlib.pyplot as plt
from matplotlib import rcParams
rcParams.update({
"font.family": "serif",
"font.serif": ["Computer Modern Roman"],
"mathtext.fontset": "cm",
"text.usetex": True, # requires LaTeX installation
"axes.linewidth": 1.2,
"axes.labelsize": 18,
"axes.titlesize": 18,
"xtick.labelsize": 16,
"ytick.labelsize": 16,
"legend.fontsize": 14,
"figure.dpi": 300,
"xtick.direction": "in",
"ytick.direction": "in",
"xtick.minor.visible": True,
"ytick.minor.visible": True,
"xtick.major.size": 6,
"xtick.minor.size": 3,
"ytick.major.size": 6,
"ytick.minor.size": 3,
"xtick.major.width": 1.0,
"ytick.major.width": 1.0,
})
In [4]:
# colors = ["tab:blue", "tab:green", "tab:orange", "tab:red", "tab:purple", "black"]
# colors = plt.cm.viridis(np.linspace(0, 1, 6))
colors = ['#EE7733', '#0077BB', '#33BBEE', '#EE3377', '#CC3311', '#009988', '#BBBBBB']
# colors = ['#D46362','#BD353A', '#5E97D0', '#256FB0', '#F2B134', '#B27D2C', '#4EB265', '#2E8540']
In [15]:
colors[::2]
Out[15]:
['#EE7733', '#33BBEE', '#CC3311', '#BBBBBB']
In [ ]:
# Draw plot concentration vs time for the old model with loops
# klin = (165,0.044,0.27,550,5000,78,4.4,5.1)
# knonlin = (0.3,2)
klin = (150, 0.04, 0.3, 500, 5000, 80, 4, 5)
knonlin = (0.3, 2)
# time_points = np.linspace(0, 11, 111)
time_points = np.linspace(0, 13, 13001) # delta t = 0.001
# initial_conditions = [10**2, 10**(-2), 10**(-2), 10**(-1)]# s_star[0]
initial_conditions = s_star[0]
# initial_conditions = [1,2,3,4]# s_star[0]
sol = solve_POSm4(initial_conditions, time_points,
klin=klin, knonlin=knonlin)
start_index = int(1.2 / 0.001)
time_ranges = time_points[start_index:]-time_points[start_index]
data_ranges = sol[start_index:, :]
plt.figure(figsize=(8, 5))
for i, c in enumerate(colors[:4]):
# plt.plot(time_points[50:], sol[50:, i], lw=2.2, color=c, label=f"S{i+1}")
plt.plot(time_ranges, data_ranges[:, i], lw=2.2, color=c, label=f"S{i+1}")
# plt.plot(time_points, sol, label=['S1', 'S2', 'S3', 'S4'], color=colors[:4])
plt.axvline(x=perturbation_point, color='gray', linestyle='--', lw=1)
plt.xlabel('Time')
plt.xlim(0, 11)
plt.xticks(np.arange(0, 12, 1))
plt.ylabel('Concentration')
plt.gca().yaxis.set_ticks_position('both')
# plt.gca().yaxis.set_tick_params(direction='inout')
# Log scale for y-axis
plt.yscale('log')
plt.title('Figure 2D: Concentration vs Time for Old Model with Loops (fixed)')
plt.legend(loc='upper left', bbox_to_anchor=(1, 1))
plt.show()
In [12]:
# Create 2x2 plots to try different parameters for the fixed model
fig, axs = plt.subplots(2, 3, figsize=(18, 10))
param_sets = [
(165, 0.044, 0.27, 550, 5000, 78, 4.4, 5.1),
(200, 0.05, 0.3, 600, 6000, 80, 5.0, 6.0),
(150, 0.04, 0.25, 500, 4000, 75, 4.0, 4.5),
(180, 0.045, 0.28, 580, 5500, 79, 4.5, 5.5),
(150, 0.04, 0.3, 500, 5000, 80, 4, 5),
(190, 0.05, 0.32, 620, 5200, 82, 5.2, 6.2)
]
for ax, params in zip(axs.flatten(), param_sets):
klin = params
sol = solve_POSm4(initial_conditions, time_points,
klin=klin, knonlin=knonlin)
for i, c in enumerate(colors[:4]):
ax.plot(time_points, sol[:, i], lw=2.2, color=c, label=f"S{i+1}")
# ax.plot(time_points, sol, label=['S1', 'S2', 'S3', 'S4'])
ax.set_xlabel('Time (a.u.)', weight='bold')
ax.set_xlim(0, 11)
ax.set_xticks(np.arange(0, 12, 1))
# Show tick marks in both left and right y-axes
ax.yaxis.set_ticks_position('both')
ax.set_ylabel('Concentration (a.u.)', weight='bold')
ax.set_yscale('log')
ax.set_title(f'Parameters: {params}')
# Make legend outside the plot
ax.legend(loc='upper left', bbox_to_anchor=(1, 1))
# ax.legend()
plt.suptitle('Figure 2D Variations: Concentration vs Time for Old Model with Loops (fixed)')
plt.tight_layout(rect=[0, 0.03, 1, 0.95])
#Show in 300 dpi
plt.gcf().set_dpi(300)
plt.show()
In [13]:
klin = (165,0.044,0.27,550,5000,78,4.4,5.1)
knonlin = (0.3,2)
# Define a helper function: odeint requires a function depending
# only on the variables (first argument) and the time, t, but
# func_POSm4 is independent from t and carries two more arguments
# (klin, knonlin).
def func_POSm4_help(x,t):
return(lde.func_POSm4(x,klin,knonlin))
# Solve the system using odeint at time points 0, 1, ..., 100,
# initial vector: (1,2,3,4)
sol = odeint(func_POSm4_help, y0 = (1,2,3,4),
t = np.linspace(0,100,11))
# We set the last point of the numeric solution as point of interest
s_star = sol[-1]
C:\Users\Admin\AppData\Local\Temp\ipykernel_26224\3783449936.py:11: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information. sol = odeint(func_POSm4_help, y0 = (1,2,3,4),
In [14]:
res_tab_2 = ld.find_loops_vset(lde.func_POSm4,vset=[s_star],
klin=klin, knonlin=knonlin,max_num_loops=10)
print(f"klin={klin}, knonlin={knonlin}")
res_tab_2['loop_rep'][0]
klin=(165, 0.044, 0.27, 550, 5000, 78, 4.4, 5.1), knonlin=(0.3, 2)
Out[14]:
| loop | length | sign | |
|---|---|---|---|
| 0 | (0, 1, 2, 3, 0) | 4 | -1 |
| 1 | (1, 2, 3, 1) | 3 | 1 |
| 2 | (0, 0) | 1 | -1 |
| 3 | (1, 1) | 1 | -1 |
| 4 | (2, 2) | 1 | -1 |
| 5 | (3, 3) | 1 | -1 |
In [22]:
# Pertubation
# klin = (150, 0.04, 0.3, 500, 5000, 80, 4, 5)
# knonlin = (0.3, 2)
delta_t = 0.001
perturbation_point = 8.6 # start perturbation at t=8.6 a.u.
time_index = int(perturbation_point/delta_t) + start_index
initial_cond_POSm4 = sol[time_index, :]
cap_range = 0.05 # time range to capture before and after perturbation
sol_1 = sol[time_index - int(cap_range/delta_t):time_index + 1, :] # from t=8.1 to t=9.1
time_points_1 = np.linspace(-cap_range, 0, int(cap_range/delta_t) + 1) # from -0.5 to 0
time_points_2 = np.linspace(0, cap_range, int(cap_range/delta_t) + 1) # from 0 to 0.5
initial_cond_perturbed = [initial_cond_POSm4[0], initial_cond_POSm4[1], initial_cond_POSm4[2], 0.75] # perturb S4 to 0.75
sol_2 = solve_POSm4(initial_cond_perturbed, time_points_2, klin=klin, knonlin=knonlin)
In [23]:
# Plot pertubation results (4 plots for S1, S2, S3, S4 in one figure)
fig, axs = plt.subplots(4, 1, figsize=(10, 10))
# Plot for S1
time_S1 = np.concatenate((time_points_1, time_points_2))
value_S1 = np.concatenate((sol_1[:, 0], sol_2[:, 0]))
axs[0].plot(time_S1, value_S1, lw=5, color=colors[0], label='S1')
# Black dashed line at t=0.0
axs[0].axvline(x=0.0, color='black', linestyle='--', lw=1.5, label='Perturbation')
axs[0].set_xlabel('Time')
axs[0].set_ylabel('Concentration of S1')
axs[0].set_title('Perturbation Analysis: Concentration of S1 over Time when S4 is perturbed')
axs[0].set_xlim(-cap_range, cap_range)
axs[0].yaxis.set_ticks_position('both')
# axs[0].set_xticks(np.arange(-cap_range, cap_range + delta_t, delta_t))
axs[0].legend(loc='upper right')
# Plot for S2
time_S2 = np.concatenate((time_points_1, time_points_2))
value_S2 = np.concatenate((sol_1[:, 1], sol_2[:, 1]))
axs[1].plot(time_S2, value_S2, lw=5, color=colors[1], label='S2')
# Black dashed line at t=0.0
axs[1].axvline(x=0.0, color='black', linestyle='--', lw=1.5, label='Perturbation')
axs[1].set_xlabel('Time')
axs[1].set_ylabel('Concentration of S2')
axs[1].set_title('Perturbation Analysis: Concentration of S2 over Time when S4 is perturbed')
axs[1].set_xlim(-cap_range, cap_range)
axs[1].yaxis.set_ticks_position('both')
# axs[1].set_xticks(np.arange(-cap_range, cap_range + delta_t, delta_t))
axs[1].legend(loc='upper right')
# Plot for S3
time_S3 = np.concatenate((time_points_1, time_points_2))
value_S3 = np.concatenate((sol_1[:, 2], sol_2[:, 2]))
axs[2].plot(time_S3, value_S3, lw=5, color=colors[2], label='S3')
# Black dashed line at t=0.0
axs[2].axvline(x=0.0, color='black', linestyle='--', lw=1.5, label='Perturbation')
axs[2].set_xlabel('Time')
axs[2].set_ylabel('Concentration of S3')
axs[2].set_title('Perturbation Analysis: Concentration of S3 over Time when S4 is perturbed')
axs[2].set_xlim(-cap_range, cap_range)
axs[2].yaxis.set_ticks_position('both')
# axs[2].set_xticks(np.arange(-cap_range, cap_range + delta_t, delta_t))
axs[2].legend(loc='upper right')
# Plot for S4
time_S4 = np.concatenate((time_points_1, time_points_2))
value_S4 = np.concatenate((sol_1[:, 3], sol_2[:, 3]))
axs[3].plot(time_S4, value_S4, lw=5, color=colors[3], label='S4')
# Black dashed line at t=0.0
axs[3].axvline(x=0.0, color='black', linestyle='--', lw=1.5, label='Perturbation')
axs[3].set_xlabel('Time')
axs[3].set_ylabel('Concentration of S4')
axs[3].set_title('Perturbation Analysis: Concentration of S4 over Time when S4 is perturbed')
axs[3].set_xlim(-cap_range, cap_range)
axs[3].yaxis.set_ticks_position('both')
# axs[3].set_xticks(np.arange(-cap_range, cap_range + delta_t, delta_t))
axs[3].legend(loc='upper right')
plt.tight_layout()
plt.show()
In [23]:
# Try solve_ivp instead of odeint
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
# Draw plot concentration vs time for the old model with loops
def model_ivp(t, x, klin, knonlin):
return lde.func_POSm4(x, klin, knonlin)
time_points = np.linspace(0, 10, 100)
initial_conditions = [10**2, 10**(-2), 10**(-2), 10**(-1)]# s_star[0]
sol_ivp = solve_ivp(model_ivp, [time_points[0], time_points[-1]], initial_conditions,
t_eval=time_points, args=(klin, knonlin), method='RK45', rtol=1e-6)
plt.plot(sol_ivp.t, sol_ivp.y.T)
plt.xlabel('Time')
plt.ylabel('Concentration')
# Log scale for y-axis
plt.yscale('log')
plt.title('Concentration vs Time for Old Model with Loops')
plt.show()
In [91]:
from matplotlib.patches import FancyArrowPatch, Rectangle
# ---------------------------------------------------------
# Helper function: draw a state box with label
# ---------------------------------------------------------
def draw_state(ax, center, label, color):
x, y = center
w, h = 0.8, 0.5
rect = Rectangle((x - w/2, y - h/2), w, h,
linewidth=1.5, edgecolor=None,
facecolor=color, alpha=0.3)
rect_2 = Rectangle((x - w/2, y - h/2), w, h,
linewidth=1.5, edgecolor="black",
fill=False)
ax.add_patch(rect)
ax.add_patch(rect_2)
ax.text(x, y, label, ha="center", va="center", fontsize=12)
return (x, y, w, h)
# ---------------------------------------------------------
# Helper: draw straight arrow between two centers
# ---------------------------------------------------------
def arrow(ax, start, end, text=None, y_offset=0.0):
x1, y1 = start
x2, y2 = end
arr = FancyArrowPatch(
(x1, y1 + y_offset), (x2, y2 + y_offset),
arrowstyle='->', mutation_scale=12, linewidth=1.5, color='black'
)
ax.add_patch(arr)
if text:
ax.text((x1+x2)/2, y1 + y_offset + 0.15, text,
ha="center", va="bottom", fontsize=12)
# ---------------------------------------------------------
# Build figure
# ---------------------------------------------------------
fig, ax = plt.subplots(figsize=(8, 3))
# Node positions (all on same horizontal line)
y0 = 0
S1 = draw_state(ax, (0, y0), "S$_1$", colors[0])
S2 = draw_state(ax, (2, y0), "S$_2$", colors[1])
S3 = draw_state(ax, (4, y0), "S$_3$", colors[2])
S4 = draw_state(ax, (6, y0), "S$_4$", colors[3])
# ---------------------------------------------------------
# Straight horizontal arrows (1→S1, S1→S2, S2→S3, S3→S4)
# ---------------------------------------------------------
arrow(ax, (-1.0, y0), (S1[0]-0.4, y0), "1") # left input
arrow(ax, (S1[0]+0.4, y0), (S2[0]-0.4, y0), "2")
arrow(ax, (S2[0]+0.4, y0), (S3[0]-0.4, y0), "4")
arrow(ax, (S3[0]+0.4, y0), (S4[0]-0.4, y0), "6")
# Vertical output arrows (3, 5, 7, 8)
for (cx, cy, _, _), label in zip([S1, S2, S3, S4], ["3", "5", "7", "8"]):
arr = FancyArrowPatch((cx, cy - 0.3), (cx, cy - 1.0),
arrowstyle='->', mutation_scale=12,
linewidth=1.5, color="black")
ax.add_patch(arr)
ax.text(cx + 0.15, cy - 0.65, label, fontsize=12, va="center")
# ---------------------------------------------------------
# Feedback arc S4 → S2 (curved)
# ---------------------------------------------------------
curve = FancyArrowPatch(
(S4[0]+0.4, y0), (S2[0]-0.4, y0),
arrowstyle='->',
mutation_scale=12,
linewidth=1.5,
connectionstyle="arc3,rad=0.5",
color="gray"
)
ax.add_patch(curve)
# Label "??" for the feedback arc (number ? )
ax.text((S4[0]+S2[0])/2, y0 + 1.0, "", fontsize=12)
# ---------------------------------------------------------
# Formatting
# ---------------------------------------------------------
ax.set_xlim(-1, 7)
ax.set_ylim(-2, 2)
ax.axis("off")
plt.tight_layout()
plt.show()
In [96]:
fig, ax = plt.subplots(figsize=(6, 3))
# Node positions (all on same horizontal line)
y0 = 0
S1 = draw_state(ax, (0, y0), "S$_1$", colors[0])
S2 = draw_state(ax, (2.5, y0), "S$_2$", colors[1])
S3 = draw_state(ax, (5.0, y0), "S$_3$", colors[2])
S4 = draw_state(ax, (7.5, y0), "S$_4$", colors[3])
# ---------------------------------------------------------
# Straight horizontal arrows (1→S1, S1→S2, S2→S3, S3→S4)
# ---------------------------------------------------------
arrow(ax, (-1.0, y0), (S1[0]-0.4, y0), "1") # left input
arrow(ax, (S1[0]+0.4, y0), (S2[0]-0.4, y0), "2")
arrow(ax, (S2[0]+0.4, y0), (S3[0]-0.4, y0), "4")
arrow(ax, (S3[0]+0.4, y0), (S4[0]-0.4, y0), "6")
# Vertical output arrows (3, 5, 7, 8)
for (cx, cy, _, _), label in zip([S1, S2, S3, S4], ["3", "5", "7", "8"]):
arr = FancyArrowPatch((cx, cy - 0.3), (cx, cy - 1.0),
arrowstyle='->', mutation_scale=12,
linewidth=1.5, color="black")
ax.add_patch(arr)
ax.text(cx + 0.15, cy - 0.65, label, fontsize=12, va="center")
# Example positions (use your own)
S1_pos = (0, 0)
S2_pos = (2.5, 0)
S3_pos = (5.0, 0)
S4_pos = (7.5, 0)
# Draw rectangular-corner feedback: S4 → S2
x4, y4 = S4_pos
x2, y2 = S2_pos
# Segment 1: Up
ax.add_patch(FancyArrowPatch(
(x4, y4+0.3), (x4, y4 + 1.03),
arrowstyle='-', linewidth=1.8, color="gray"))
# Segment 2: Left
ax.add_patch(FancyArrowPatch(
(x4+0.05, y4 + 1), (x2-0.05, y2 + 1),
arrowstyle='-', linewidth=1.8, color="gray"))
# Segment 3: Down with arrow head
ax.add_patch(FancyArrowPatch(
(x2, y2 + 1.03), (x2, y2+0.25),
arrowstyle='->', mutation_scale=14, linewidth=1.8, color="gray"))
# ax.set_xlim(-1, 10)
# ax.set_ylim(-2, 3)
# ax.axis("off")
# plt.show()
# Label "??" for the feedback arc (number ? )
ax.text((S4[0]+S2[0])/2, y0 + 1.0, "", fontsize=12)
# ---------------------------------------------------------
# Formatting
# ---------------------------------------------------------
ax.set_xlim(-2, 9)
ax.set_ylim(-2, 2)
ax.axis("off")
plt.tight_layout()
plt.show()
Example 2: Model of complex formation¶
The model of complex formation has also been investigated in [4]. It is a basic model of a bilinear reaction in that two species, A, B, reversibly form a complex AB together. Only the two reactions of complex formation and complex decomposition are modelled. The model contains three species, two rate coefficients and no other parameters.
$$\frac{dA}{dt} = -k_1\,A\,B + k_2\,AB$$
$$\frac{dB}{dt} = -k_1\,A\,B + k_2\,AB$$
$$\frac{d(AB)}{dt} = k_1\,A\,B - k_2\,AB$$
We employed the following parameter set: $k_1 = 1$, $k_2 = 2$, with initial conditions $A = 10, B = 5, C = 2$ (C should change to AB).
In [38]:
def model_complex_formation(y, params):
A, B, AB = y
k1, k2 = params
dAdt = -k1 * A * B + k2 * AB
dBdt = -k1 * A * B + k2 * AB
dABdt = k1 * A * B - k2 * AB
return np.array([dAdt, dBdt, dABdt], dtype=float)
def solve_complex_formation(initial_conditions, time_points, k1, k2):
def model(t, x):
return model_complex_formation(x, (k1, k2))
sol = solve_ivp(model, [time_points[0], time_points[-1]], initial_conditions,
t_eval=time_points, method='RK45', rtol=1e-6)
return sol
In [39]:
# Detect loops in the complex formation model
s_star = [[1.0, 1.0, 1.0]]
res_tab_complex = ld.find_loops_vset(model_complex_formation,
vset=s_star,
params=[1.0, 2.0],
max_num_loops=10)
print(res_tab_complex['loop_rep'][0])
loop length sign 0 (0, 1, 0) 2 1 1 (0, 1, 2, 0) 3 -1 2 (0, 2, 0) 2 1 3 (0, 2, 1, 0) 3 -1 4 (1, 2, 1) 2 1 5 (0, 0) 1 -1 6 (1, 1) 1 -1 7 (2, 2) 1 -1
In [40]:
pd.DataFrame(res_tab_complex['loop_rep'][0]).sort_values(by=['length','loop'])
Out[40]:
| loop | length | sign | |
|---|---|---|---|
| 5 | (0, 0) | 1 | -1 |
| 6 | (1, 1) | 1 | -1 |
| 7 | (2, 2) | 1 | -1 |
| 0 | (0, 1, 0) | 2 | 1 |
| 2 | (0, 2, 0) | 2 | 1 |
| 4 | (1, 2, 1) | 2 | 1 |
| 1 | (0, 1, 2, 0) | 3 | -1 |
| 3 | (0, 2, 1, 0) | 3 | -1 |
In [41]:
id_to_subject = {
0: "A",
1: "B",
2: "AB"
}
df = pd.DataFrame(res_tab_complex['loop_rep'][0]).sort_values(by=['length','loop'])
df["loop_subject"] = df["loop"].apply(
lambda t: tuple(id_to_subject[i] for i in t)
)
df
Out[41]:
| loop | length | sign | loop_subject | |
|---|---|---|---|---|
| 5 | (0, 0) | 1 | -1 | (A, A) |
| 6 | (1, 1) | 1 | -1 | (B, B) |
| 7 | (2, 2) | 1 | -1 | (AB, AB) |
| 0 | (0, 1, 0) | 2 | 1 | (A, B, A) |
| 2 | (0, 2, 0) | 2 | 1 | (A, AB, A) |
| 4 | (1, 2, 1) | 2 | 1 | (B, AB, B) |
| 1 | (0, 1, 2, 0) | 3 | -1 | (A, B, AB, A) |
| 3 | (0, 2, 1, 0) | 3 | -1 | (A, AB, B, A) |
In [42]:
# Calculate Jacobian at s_star
j_matrix = numdifftools.Jacobian(model_complex_formation, method='central')(np.array(s_star[0]), params=[1.0, 2.0])
signed_jacobian = np.sign(j_matrix)
print("Signed Jacobian matrix at s_star:")
print(signed_jacobian)
Signed Jacobian matrix at s_star: [[-1. -1. 1.] [-1. -1. 1.] [ 1. 1. -1.]]
In [ ]:
# Visualize Jacobian Matrix as correlation matrix (Matplotlib/ Seaborn/ ...) --> Can see both name tag (which vs which) + pos/neg relationship
In [30]:
time_points = np.linspace(0, 1, 1111) # delta t = 0.01
initial_conditions = [10.0, 5.0, 2.0] # Initial concentrations of A, B, and AB
k1 = 1.0 # Rate constant for complex formation
k2 = 2.0 # Rate constant for complex dissociation
sol = solve_complex_formation(initial_conditions, time_points, k1, k2)
In [31]:
plt.figure(figsize=(8, 6))
for i, c in enumerate(colors[:3]):
plt.plot(time_points, sol.y[i], lw=5, color=c, label=f"S{i+1}")
# plt.plot(sol.t, sol.y.T)
plt.axvline(x=0.2, color='black', linestyle='--', lw=1.5, alpha=0.7)
plt.xlim(0, 1)
plt.xlabel('Time')
plt.ylabel('Concentration')
plt.title('Concentration vs Time for Complex Formation Model (Ex. 2)')
plt.legend(['A', 'B', 'AB'])
plt.show()
In [77]:
# Take initial conditions 1 at t = 0.2 of the the above simulation
initial_conditions_1 = sol.y[:, 20] # Corresponds to t = 0.2
In [43]:
# Pertubation analysis for complex formation model (Perturb A - fixed)
delta_t = 0.001
perturbation_point = 0.2 # start perturbation at t=0.2
time_index = int(perturbation_point/delta_t)
cap_range = 0.2 # time range to capture before and after perturbation
time_points_1 = np.linspace(-cap_range, 0, int(cap_range/delta_t) + 1) # delta 1= 0.001
initial_conditions_1 = sol.y[:, time_index] # Initial concentrations of A, B, and AB
# Take partial values for sol_1 at the previous simulation
sol_1_005 = sol.y[:, time_index - int(cap_range/delta_t):time_index + 1] # Corresponds to t = -0.05 to t = 0.0
time_points_2 = np.linspace(0, cap_range, int(cap_range/delta_t) + 1) # delta 2= 0.001
initial_conditions_2 = [40.0, sol.y[1, time_index], sol.y[2, time_index]] # Initial concentrations of A, B, and AB
sol_2 = solve_complex_formation(initial_conditions_2, time_points_2, k1, k2)
In [44]:
# Plot pertubation results (3 plots for A, B, and C in one figure)
fig, axs = plt.subplots(3, 1, figsize=(10, 8))
# Plot for A
time_A = np.concatenate((time_points_1, time_points_2))
value_A = np.concatenate((sol_1_005[0], sol_2.y[0]))
axs[0].plot(time_A, value_A, lw=5, color=colors[0], label='A (to 40)')
# Black dashed line at t=0.0
axs[0].axvline(x=0.0, color='black', linestyle='--', lw=1.5, label='Perturbation')
axs[0].set_xlabel('Time')
axs[0].set_ylabel('Concentration of A')
axs[0].set_title('Perturbation Analysis: Concentration of A over Time when A is perturbed')
axs[0].set_xlim(-cap_range, cap_range)
axs[0].yaxis.set_ticks_position('both')
# axs[0].set_xticks(np.arange(-cap_range, cap_range, 0.01))
axs[0].legend(loc='lower right')
# Plot for B
time_B = np.concatenate((time_points_1, time_points_2))
value_B = np.concatenate((sol_1_005[1], sol_2.y[1]))
axs[1].plot(time_B, value_B, lw=5, color=colors[1], label='B') # (from 5.0 to 5.0)
# Black dashed line at t=0.0
axs[1].axvline(x=0.0, color='black', linestyle='--', lw=1.5, label='Perturbation')
axs[1].set_xlabel('Time')
axs[1].set_ylabel('Concentration of B')
axs[1].set_title('Perturbation Analysis: Concentration of B over Time when A is perturbed')
axs[1].set_xlim(-cap_range, cap_range)
axs[1].yaxis.set_ticks_position('both')
# axs[1].set_xticks(np.arange(-cap_range, cap_range, 0.01))
axs[1].legend(loc='upper right')
# Plot for AB
time_AB = np.concatenate((time_points_1, time_points_2))
value_AB = np.concatenate((sol_1_005[2], sol_2.y[2]))
axs[2].plot(time_AB, value_AB, lw=5, color=colors[2], label='AB') # (from 2.0 to 2.0)
# Black dashed line at t=0.0
axs[2].axvline(x=0.0, color='black', linestyle='--', lw=1.5, label='Perturbation')
axs[2].set_xlabel('Time')
axs[2].set_ylabel('Concentration of AB')
axs[2].set_title('Perturbation Analysis: Concentration of AB over Time when A is perturbed')
axs[2].set_xlim(-cap_range, cap_range)
axs[2].yaxis.set_ticks_position('both')
# axs[2].set_xticks(np.arange(-cap_range, cap_range, 0.01))
axs[2].legend(loc='lower right')
plt.tight_layout()
plt.show()
In [100]:
# Pertubation analysis for complex formation model (Perturb B)
time_points_1 = np.linspace(-0.05, 0, 51) # delta 1= 0.001
initial_conditions_1 = sol.y[:, 200] # Initial concentrations of A, B, and AB
# Take partial values for sol_1 at the previous simulation
sol_1_005 = sol.y[:, 150:201] # Corresponds to t = -0.05 to t = 0.0
time_points_2 = np.linspace(0, 0.05, 51) # delta 2= 0.001
initial_conditions_2 = [sol.y[0, 200], 10.0, sol.y[2, 200]] # Initial concentrations of A, B, and AB
sol_2 = solve_complex_formation(initial_conditions_2, time_points_2, k1, k2)
In [102]:
# Plot pertubation results (3 plots for A, B, and C in one figure)
fig, axs = plt.subplots(3, 1, figsize=(10, 8))
# Plot for A
time_A = np.concatenate((time_points_1, time_points_2))
value_A = np.concatenate((sol_1_005[0], sol_2.y[0]))
axs[0].plot(time_A, value_A, lw=5, color=colors[0], label='A')
# Black dashed line at t=0.0
axs[0].axvline(x=0.0, color='black', linestyle='--', lw=1.5, label='Perturbation')
axs[0].set_xlabel('Time')
axs[0].set_ylabel('Concentration of A')
axs[0].set_title('Perturbation Analysis: Concentration of A over Time when B is perturbed')
axs[0].set_xlim(-0.05, 0.05)
axs[0].yaxis.set_ticks_position('both')
axs[0].set_xticks(np.arange(-0.05, 0.05, 0.01))
axs[0].legend(loc='upper right')
# Plot for B
time_B = np.concatenate((time_points_1, time_points_2))
value_B = np.concatenate((sol_1_005[1], sol_2.y[1]))
axs[1].plot(time_B, value_B, lw=5, color=colors[1], label='B (to 10.0)') # (from 5.0 to 5.0)
# Black dashed line at t=0.0
axs[1].axvline(x=0.0, color='black', linestyle='--', lw=1.5, label='Perturbation')
axs[1].set_xlabel('Time')
axs[1].set_ylabel('Concentration of B')
axs[1].set_title('Perturbation Analysis: Concentration of B over Time when B is perturbed')
axs[1].set_xlim(-0.05, 0.05)
axs[1].yaxis.set_ticks_position('both')
axs[1].set_xticks(np.arange(-0.05, 0.05, 0.01))
axs[1].legend(loc='lower right')
# Plot for AB
time_AB = np.concatenate((time_points_1, time_points_2))
value_AB = np.concatenate((sol_1_005[2], sol_2.y[2]))
axs[2].plot(time_AB, value_AB, lw=5, color=colors[2], label='AB') # (from 2.0 to 2.0)
# Black dashed line at t=0.0
axs[2].axvline(x=0.0, color='black', linestyle='--', lw=1.5, label='Perturbation')
axs[2].set_xlabel('Time')
axs[2].set_ylabel('Concentration of AB')
axs[2].set_title('Perturbation Analysis: Concentration of AB over Time when B is perturbed')
axs[2].set_xlim(-0.05, 0.05)
axs[2].yaxis.set_ticks_position('both')
axs[2].set_xticks(np.arange(-0.05, 0.05, 0.01))
axs[2].legend(loc='lower right')
plt.tight_layout()
plt.show()
Example 3: Model of complex formation with activation-deactivation cycle¶
The model of complex formation with activation-deactivation cycle has also been investigated in [4]. It is composed of the basic model of complex formation in that A and B bind to form complex AB as described in section 1.2, together with an activation and deactivation of one of the complex-forming species, B and B∗. In addition to the two reactions of complex formation and complex decomposition, the activating reaction and the deactivation reaction are modelled using Michaelis-Menten reaction kinetics. Thereby, the activation of B to B∗ is in addition stimulus-dependent (the stimulus is modelled as a parameter, kn3), while the deactivation of B∗ to B depends linearly on species A. The model contains four species, A, B, AB, B∗ , four rate coefficients, k1 . . . , k4, two activation constants, kn1, kn2, and a stimulus parameter, kn3.
$$ \begin{aligned} \frac{dA}{dt} &= -k_1\,A\,B + k_2\,AB ,\\[6pt] \frac{dB}{dt} &= -k_1\,A\,B + k_2\,AB + k_3\,A\,\frac{B^*}{B^* + kn_1} - k_4\,kn_3\,\frac{B}{B + kn_2} ,\\[6pt] \frac{d(AB)}{dt} &= k_1\,A\,B - k_2\,AB ,\\[6pt] \frac{dB^*}{dt} &= -k_3\,A\,\frac{B^*}{B^* + kn_1} + k_4\,kn_3\,\frac{B}{B + kn_2} . \end{aligned} $$
We employed the following parameter set for the analyses: k1 = 0.021458, k2 = 0.000227335, k3 = 663.297, k4 = 0.000139797, kn1 = 55.5913, kn2 = 10.2194, and values for the stimulus of kn3 = 1 (unstimulated) or kn3 = 20 (stimulated). The initial conditions were set as A = 5, B = 300, AB = 0.0001, B∗ = 0.01.
In [18]:
def model_complex_formation_act_deact_cycle(y, params):
# A, B, AB, B_star = y
# print(y[0])
A = y[0]
B = y[1]
AB = y[2]
B_star = y[3]
# print(f"params: {params}")
k1, k2, k3, k4, kn1, kn2, kn3 = params
dAdt, dBdt, dABdt, dB_stardt = 0,0,0,0 # I'm curious if it is redundant to do this step
dAdt = -k1 * A * B + k2 * AB
dBdt = -k1 * A * B + k2 * AB + k3 * A * (B_star / (B_star + kn1)) - k4 * kn3 * (B / (B + kn2))
dABdt = k1 * A * B - k2 * AB
dB_stardt = -k3 * A * (B_star / (B_star + kn1)) + k4 * kn3 * (B / (B + kn2))
# return [dAdt, dBdt, dABdt, dB_stardt]
return np.array([dAdt, dBdt, dABdt, dB_stardt], dtype=float)
def solve_complex_formation_act_deact_cycle(initial_conditions, time_points, params):
def model(t, x):
return model_complex_formation_act_deact_cycle(x, params)
sol = solve_ivp(model, [time_points[0], time_points[-1]], initial_conditions,
t_eval=time_points, method='RK45', rtol=1e-6)
return sol
In [19]:
# Initial conditions and parameters for the activation-deactivation cycle model
initial_conditions_act_deact = [10.0, 5.0, 2.0, 3.0] # Initial concentrations of A, B, AB, and B*
k1 = 1.0 # Rate constant for complex formation
k2 = 2.0 # Rate constant for complex dissociation
k3 = 0.5 # Rate constant for activation
k4 = 0.3 # Rate constant for deactivation
kn1 = 1.0 # Michaelis constant for activation
kn2 = 1.0 # Michaelis constant for deactivation
kn3 = 1.0 # Scaling factor for deactivation
params_act_deact = [k1, k2, k3, k4, kn1, kn2, kn3]
time_points_act_deact = np.linspace(0, 1, 111)
In [20]:
# Solve the activation-deactivation cycle model
sol_act_deact = solve_complex_formation_act_deact_cycle(initial_conditions_act_deact,
time_points_act_deact,
[k1, k2, k3, k4, kn1, kn2, kn3])
In [21]:
plt.figure(figsize=(8, 6))
for i, c in enumerate(colors[:4]):
plt.plot(time_points_act_deact, sol_act_deact.y[i], lw=5, color=c, label=f"S{i+1}")
# plt.plot(sol_act_deact.t, sol_act_deact.y.T)
plt.xlim(0, 1)
plt.xlabel('Time')
plt.ylabel('Concentration')
plt.title('Concentration vs Time for Activation-Deactivation Cycle Model')
plt.legend(['A', 'B', 'AB', 'B*'])
plt.show()
In [22]:
# Loop detection for activation-deactivation cycle model
s_star_3 = [(1.0, 1.0, 1.0, 1.0)]
res_tab_act_deact = ld.find_loops_vset(model_complex_formation_act_deact_cycle,
vset=s_star_3,
params=params_act_deact,
max_num_loops=10)
print(res_tab_act_deact['loop_rep'][0])
loop length sign 0 (0, 1, 0) 2 1 1 (0, 1, 2, 0) 3 -1 2 (0, 2, 0) 2 1 3 (0, 2, 1, 0) 3 -1 4 (0, 3, 1, 0) 3 1 5 (0, 3, 1, 2, 0) 4 -1 6 (1, 2, 1) 2 1 7 (1, 3, 1) 2 1 8 (0, 0) 1 -1 9 (1, 1) 1 -1
In [57]:
ld.sort_loop_index(res_tab_act_deact['loop_rep'][0])
Out[57]:
| loop | length | sign | |
|---|---|---|---|
| 0 | (0, 1, 0) | 2 | 1 |
| 1 | (0, 1, 2, 0) | 3 | -1 |
| 2 | (0, 2, 0) | 2 | 1 |
| 3 | (0, 2, 1, 0) | 3 | -1 |
| 4 | (0, 3, 1, 0) | 3 | 1 |
| 5 | (0, 3, 1, 2, 0) | 4 | -1 |
| 6 | (1, 2, 1) | 2 | 1 |
| 7 | (1, 3, 1) | 2 | 1 |
| 8 | (0, 0) | 1 | -1 |
| 9 | (1, 1) | 1 | -1 |
In [24]:
id_to_subject_act_deact = {
0: "A",
1: "B",
2: "AB",
3: "B*"
}
df_act_deact = pd.DataFrame(res_tab_act_deact['loop_rep'][0]).sort_values(by=['length','loop'])
df_act_deact["loop_subject"] = df_act_deact["loop"].apply(
lambda t: tuple(id_to_subject_act_deact[i] for i in t)
)
df_act_deact
Out[24]:
| loop | length | sign | loop_subject | |
|---|---|---|---|---|
| 8 | (0, 0) | 1 | -1 | (A, A) |
| 9 | (1, 1) | 1 | -1 | (B, B) |
| 0 | (0, 1, 0) | 2 | 1 | (A, B, A) |
| 2 | (0, 2, 0) | 2 | 1 | (A, AB, A) |
| 6 | (1, 2, 1) | 2 | 1 | (B, AB, B) |
| 7 | (1, 3, 1) | 2 | 1 | (B, B*, B) |
| 1 | (0, 1, 2, 0) | 3 | -1 | (A, B, AB, A) |
| 3 | (0, 2, 1, 0) | 3 | -1 | (A, AB, B, A) |
| 4 | (0, 3, 1, 0) | 3 | 1 | (A, B*, B, A) |
| 5 | (0, 3, 1, 2, 0) | 4 | -1 | (A, B*, B, AB, A) |
In [27]:
# Calculate Jacobian matrix at s_star_3
j_matrix = numdifftools.Jacobian(
lambda y: model_complex_formation_act_deact_cycle(y, params_act_deact),
method="central"
)([1,1,1,1])
signed_jacobian = np.sign(j_matrix)
In [28]:
signed_jacobian
Out[28]:
array([[-1., -1., 1., 0.],
[-1., -1., 1., 1.],
[ 1., 1., -1., 0.],
[-1., 1., 0., -1.]])
In [29]:
loop_list = ld.find_loops(j_matrix)
ld.loop_summary(loop_list)
Out[29]:
| length | 1 | 2 | 3 | 4 |
|---|---|---|---|---|
| total | 4 | 4 | 3 | 1 |
| neg | 4 | 0 | 2 | 1 |
| pos | 0 | 4 | 1 | 0 |
In [ ]:
# Loop summary (A different way)
loop_summary = df_act_deact.groupby('length').agg(
num_loops=('loop', 'count'),
loops=('loop_subject', lambda x: list(x))
).reset_index()
loop_summary
Out[ ]:
| length | num_loops | loops | |
|---|---|---|---|
| 0 | 1 | 2 | [(A, A), (B, B)] |
| 1 | 2 | 4 | [(A, B, A), (A, AB, A), (B, AB, B), (B, B*, B)] |
| 2 | 3 | 3 | [(A, B, AB, A), (A, AB, B, A), (A, B*, B, A)] |
| 3 | 4 | 1 | [(A, B*, B, AB, A)] |
Example 4: Full MAPK cascade model¶
In [ ]:
def model_MAPK_22(y, k):
"""
22-variable MAPK enzymatic model.
Variables: S1 ... S22
Parameters: k1 ... k30 passed as a single list or array k[0] ... k[29]
"""
# Unpack variables
S1, S2, S3, S4, S5, S6, S7, S8, S9, S10, \
S11, S12, S13, S14, S15, S16, S17, S18, S19, S20, \
S21, S22 = y
# Unpack parameters (k1 ... k30)
k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, \
k11, k12, k13, k14, k15, k16, k17, k18, k19, k20, \
k21, k22, k23, k24, k25, k26, k27, k28, k29, k30 = k
# ODE system (exactly as in the paper)
dS1 = -k1*S1*S3 + k2*S13 + k3*S13
dS2 = -k4*S2*S4 + k5*S14 + k6*S14
dS3 = -k1*S1*S3 + k2*S13 + k6*S14
dS4 = k3*S13 - k4*S2*S4 + k5*S14 \
- k7*S5*S4 + k8*S15 + k9*S15
dS5 = -k7*S5*S4 + k8*S15 + k12*S20
dS6 = k9*S15 - k10*S6*S12 + k11*S20 \
- k13*S6*S4 + k14*S16 + k18*S19
dS7 = k15*S16 - k10*S7*k17 + k17*S19 \
- k19*S8*S7 + k20*S18 + k27*S18
dS8 = -k19*S8*S7 + k20*S17 + k24*S22
dS9 = k21*S17 - k22*S9*S11 + k23*S22 \
- k25*S9*S7 + k26*S18 + k30*S21
dS10 = k27*S18 - k28*S10*S11 + k29*S21
dS11 = -k22*S9*S11 + k23*S22 - k28*S10*S11 + k29*S21
dS12 = -k10*S6*S12 + k11*S20 - k17*S12 + k18*S19
dS13 = k1*S1*S3 - k2*S13 - k3*S13
dS14 = k4*S2*S4 - k5*S14 - k6*S14
dS15 = k7*S5*S4 - k8*S15 - k9*S15
dS16 = k13*S6*S4 - k14*S16 - k15*S16
dS17 = k19*S8*S7 - k20*S17 - k21*S17
dS18 = k25*S9*S7 - k26*S18 - k27*S18
dS19 = k16*S7*S12 - k17*S19 - k18*S19
dS20 = k10*S6*S12 - k11*S20 - k12*S20
dS21 = k28*S10*S11 - k29*S21 - k30*S21
dS22 = k22*S9*S11 - k23*S22 - k24*S22
return np.array([
dS1, dS2, dS3, dS4, dS5, dS6, dS7, dS8, dS9, dS10,
dS11, dS12, dS13, dS14, dS15, dS16, dS17, dS18,
dS19, dS20, dS21, dS22
], dtype=float)
def solve_MAPK_22(initial_conditions, time_points, k):
def ode(t, x):
return model_MAPK_22(x, k)
sol = solve_ivp(
ode,
[time_points[0], time_points[-1]],
initial_conditions,
t_eval=time_points,
method="RK45",
rtol=1e-6
)
return sol
In [ ]:
In [61]:
res_tab = ld.find_loops_vset(lde.func_li08,
vset=np.ones((1,18)),
t=0,
max_num_loops=100)
# The loop list, a pandas dataframe, is accessed like this:
res_tab['loop_rep'][0]
Out[61]:
| loop | length | sign | |
|---|---|---|---|
| 0 | (0, 1, 0) | 2 | -1 |
| 1 | (0, 2, 1, 0) | 3 | 1 |
| 2 | (0, 3, 4, 6, 0) | 4 | -1 |
| 3 | (0, 8, 9, 10, 2, 1, 0) | 6 | 1 |
| 4 | (0, 8, 9, 12, 0) | 4 | -1 |
| 5 | (0, 8, 9, 13, 3, 4, 6, 0) | 7 | 1 |
| 6 | (8, 9, 11, 8) | 3 | -1 |
| 7 | (1, 2, 1) | 2 | -1 |
| 8 | (0, 0) | 1 | -1 |
| 9 | (1, 1) | 1 | -1 |
| 10 | (2, 2) | 1 | -1 |
| 11 | (3, 3) | 1 | -1 |
| 12 | (4, 4) | 1 | -1 |
| 13 | (5, 5) | 1 | -1 |
| 14 | (6, 6) | 1 | -1 |
| 15 | (7, 7) | 1 | -1 |
| 16 | (8, 8) | 1 | -1 |
| 17 | (9, 9) | 1 | -1 |
| 18 | (10, 10) | 1 | -1 |
| 19 | (11, 11) | 1 | -1 |
| 20 | (12, 12) | 1 | -1 |
| 21 | (13, 13) | 1 | -1 |
| 22 | (15, 15) | 1 | 1 |
In [59]:
# Jacobian matrix at the point of interest
j_matrix_li08 = numdifftools.Jacobian(lde.func_li08)(
np.ones((18,)),
t=0
).real
signed_jacobian_li08 = np.sign(j_matrix_li08)
print("Signed Jacobian matrix for li08 model at point of interest:")
print(signed_jacobian_li08)
Signed Jacobian matrix for li08 model at point of interest: [[-1. 1. 0. 0. 0. 0. -1. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 0.] [-1. -1. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 1. -1. -1. 0. 0. 0. 0. 0. 0. 0. -1. 0. 0. 0. 0. 0. 0. 0.] [ 1. 0. 0. -1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0.] [ 0. 0. 0. -1. -1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 1. 0. 0. 0. 1. -1. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. -1. 0. -1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 1. 0. 0. 0. 0. 0. 0. -1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 1. 0. 0. 0. 0. 0. 0. 0. -1. 0. 0. 1. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 1. -1. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. -1. -1. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. -1. 0. -1. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. -1. 0. 0. -1. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. -1. 0. 0. 0. -1. 0. 0. 0. 0.] [-1. 1. 1. 0. 0. 0. 0. 0. 0. 0. -1. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 1.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 1.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]
In [60]:
loop_list_li08 = ld.find_loops(j_matrix_li08)
ld.loop_summary(loop_list_li08)
Out[60]:
| length | 1 | 2 | 3 | 4 | 6 | 7 |
|---|---|---|---|---|---|---|
| total | 15 | 2 | 2 | 2 | 1 | 1 |
| neg | 14 | 2 | 1 | 2 | 0 | 0 |
| pos | 1 | 0 | 1 | 0 | 1 | 1 |
In [58]:
# Jacobian matrix at the point of interest
j_matrix_li08 = numdifftools.Jacobian(lde.func_li08,method="complex")(
np.ones((18,)),
t=0
).real
signed_jacobian_li08 = np.sign(j_matrix_li08)
print("Signed Jacobian matrix for li08 model at point of interest:")
print(signed_jacobian_li08)
Signed Jacobian matrix for li08 model at point of interest: [[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]
In [ ]:
# Solve MAPK cascade model (li08) and plot results
def model_li08(t, x):
return lde.func_li08(x, t)
time_points_li08 = np.linspace(0, 11, 111)
initial_conditions_li08 = np.ones((18,)) # Initial concentrations for li08 model
sol_li08 = solve_ivp(model_li08, [time_points_li08[0], time_points_li08[-1]], initial_conditions_li08,
t_eval=time_points_li08, method='RK45', rtol=1e-6)
In [109]:
plt.figure(figsize=(10, 6))
plt.plot(sol_li08.t, sol_li08.y.T)
plt.xlim(0, 11)
plt.xlabel('Time')
plt.ylabel('Concentration')
plt.title('Concentration vs Time for MAPK Cascade Model (li08)')
plt.legend([f'S{i+1}' for i in range(18)], bbox_to_anchor=(1.05, 1), loc='upper left')
plt.tight_layout()
plt.show()
Reimplementation¶
In [87]:
def model_MAPK_22(y, k):
"""
22-variable MAPK enzymatic model.
Variables: S1 ... S22
Parameters: k1 ... k30 passed as a single list or array k[0] ... k[29]
"""
# Unpack variables
S1, S2, S3, S4, S5, S6, S7, S8, S9, S10, \
S11, S12, S13, S14, S15, S16, S17, S18, S19, S20, \
S21, S22 = y
# Unpack parameters (k1 ... k30)
k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, \
k11, k12, k13, k14, k15, k16, k17, k18, k19, k20, \
k21, k22, k23, k24, k25, k26, k27, k28, k29, k30 = k
# ODE system (exactly as in the paper)
dS1 = -k1*S1*S3 + k2*S13 + k3*S13
dS2 = -k4*S2*S4 + k5*S14 + k6*S14
dS3 = -k1*S1*S3 + k2*S13 + k6*S14
dS4 = k3*S13 - k4*S2*S4 + k5*S14 \
- k7*S5*S4 + k8*S15 + k9*S15
dS5 = -k7*S5*S4 + k8*S15 + k12*S20
dS6 = k9*S15 - k10*S6*S12 + k11*S20 \
- k13*S6*S4 + k14*S16 + k18*S19
dS7 = k15*S16 - k10*S7*k17 + k17*S19 \
- k19*S8*S7 + k20*S18 + k27*S18
dS8 = -k19*S8*S7 + k20*S17 + k24*S22
dS9 = k21*S17 - k22*S9*S11 + k23*S22 \
- k25*S9*S7 + k26*S18 + k30*S21
dS10 = k27*S18 - k28*S10*S11 + k29*S21
dS11 = -k22*S9*S11 + k23*S22 - k28*S10*S11 + k29*S21
dS12 = -k10*S6*S12 + k11*S20 - k17*S12 + k18*S19
dS13 = k1*S1*S3 - k2*S13 - k3*S13
dS14 = k4*S2*S4 - k5*S14 - k6*S14
dS15 = k7*S5*S4 - k8*S15 - k9*S15
dS16 = k13*S6*S4 - k14*S16 - k15*S16
dS17 = k19*S8*S7 - k20*S17 - k21*S17
dS18 = k25*S9*S7 - k26*S18 - k27*S18
dS19 = k16*S7*S12 - k17*S19 - k18*S19
dS20 = k10*S6*S12 - k11*S20 - k12*S20
dS21 = k28*S10*S11 - k29*S21 - k30*S21
dS22 = k22*S9*S11 - k23*S22 - k24*S22
return np.array([
dS1, dS2, dS3, dS4, dS5, dS6, dS7, dS8, dS9, dS10,
dS11, dS12, dS13, dS14, dS15, dS16, dS17, dS18,
dS19, dS20, dS21, dS22
], dtype=float)
def solve_MAPK_22(initial_conditions, time_points, k):
def ode(t, x):
return model_MAPK_22(x, k)
sol = solve_ivp(
ode,
[time_points[0], time_points[-1]],
initial_conditions,
t_eval=time_points,
method="RK45",
rtol=1e-6
)
return sol
In [88]:
# Activation rates
big = 1000.0
# Dephos / dissociation rates
mid = 150.0
params_22 = [
big, mid, mid, # k1, k2, k3
big, mid, mid, # k4, k5, k6
big, mid, mid, # k7, k8, k9
big, mid, mid, # k10, k11, k12
big, mid, mid, # k13, k14, k15
big, mid, mid, # k16, k17, k18
big, mid, mid, # k19, k20, k21
big, mid, mid, # k22, k23, k24
big, mid, mid, # k25, k26, k27
big, mid, mid # k28, k29, k30
]
initial_conditions_22 = [
3e-5, # S1 E1
3e-4, # S2 E2
0.003, # S3 MAPKKK
0.0, # S4
1.2, # S5 MAPKK
0.0, # S6
0.0, # S7
1.2, # S8 MAPK
0.0, # S9
0.0, # S10
0.12, # S11 KPase
3e-4, # S12 KK_Pase
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 # S13–S22
]
In [89]:
t = np.linspace(0, 2000, 20001)
sol = solve_MAPK_22(initial_conditions_22, t, params_22)
--------------------------------------------------------------------------- KeyboardInterrupt Traceback (most recent call last) Cell In[89], line 2 1 t = np.linspace(0, 2000, 20001) ----> 2 sol = solve_MAPK_22(initial_conditions_22, t, params_22) Cell In[87], line 56, in solve_MAPK_22(initial_conditions, time_points, k) 53 def ode(t, x): 54 return model_MAPK_22(x, k) ---> 56 sol = solve_ivp( 57 ode, 58 [time_points[0], time_points[-1]], 59 initial_conditions, 60 t_eval=time_points, 61 method="RK45", 62 rtol=1e-6 63 ) 64 return sol File d:\Work\StudentJob\Jana Wolf\LoopDetect_2025\env\lib\loopdetect_python_3.12\Lib\site-packages\scipy\integrate\_ivp\ivp.py:655, in solve_ivp(fun, t_span, y0, method, t_eval, dense_output, events, vectorized, args, **options) 653 status = None 654 while status is None: --> 655 message = solver.step() 657 if solver.status == 'finished': 658 status = 0 File d:\Work\StudentJob\Jana Wolf\LoopDetect_2025\env\lib\loopdetect_python_3.12\Lib\site-packages\scipy\integrate\_ivp\base.py:197, in OdeSolver.step(self) 195 else: 196 t = self.t --> 197 success, message = self._step_impl() 199 if not success: 200 self.status = 'failed' File d:\Work\StudentJob\Jana Wolf\LoopDetect_2025\env\lib\loopdetect_python_3.12\Lib\site-packages\scipy\integrate\_ivp\rk.py:144, in RungeKutta._step_impl(self) 141 h = t_new - t 142 h_abs = np.abs(h) --> 144 y_new, f_new = rk_step(self.fun, t, y, self.f, h, self.A, 145 self.B, self.C, self.K) 146 scale = atol + np.maximum(np.abs(y), np.abs(y_new)) * rtol 147 error_norm = self._estimate_error_norm(self.K, h, scale) File d:\Work\StudentJob\Jana Wolf\LoopDetect_2025\env\lib\loopdetect_python_3.12\Lib\site-packages\scipy\integrate\_ivp\rk.py:63, in rk_step(fun, t, y, f, h, A, B, C, K) 61 K[0] = f 62 for s, (a, c) in enumerate(zip(A[1:], C[1:]), start=1): ---> 63 dy = np.dot(K[:s].T, a[:s]) * h 64 K[s] = fun(t + c * h, y + dy) 66 y_new = y + h * np.dot(K[:-1].T, B) KeyboardInterrupt:
In [ ]:
fig, ax = plt.subplots(figsize=(8, 6))
ax.plot(sol.t, sol.y[7, :], label="MAPK (S8)", color="blue", lw=2)
ax.plot(sol.t, sol.y[6, :], label="MAPK-P (S9)", color="orange", lw=2)
ax.plot(sol.t, sol.y[5, :], label="MAPK-PP (S10)", color="green", lw=2)
ax.set_xlabel("Time")
ax.set_ylabel("Concentration")
ax.set_title("Concentration vs Time for 22-Variable MAPK Model")
ax.legend()
Example 5: Reduced MAPK cascade model with and without explicit feedback regulation¶
In [68]:
def model_mapk_simple(y, klin, kn, n1):
S1, S2, S3, S4, S5, S6, S7, S8 = y
# Unpack kinetic parameters
k1, k2, k3, k4, k5, k6, k7, k8, k9, k10 = klin
# Unpack saturation constants
kn1, kn2, kn3, kn4, kn5, kn6, kn7, kn8, kn9, kn10, kn11 = kn
# dS1/dt
term_inhibit = (kn1**n1) / (kn1**n1 + S8**n1)
dS1 = -k1 * (S1/(kn2 + S1)) * term_inhibit \
+ k2 * (S2/(S2 + kn3))
# dS2/dt
dS2 = k1 * (S1/(kn2 + S1)) * term_inhibit \
- k2 * (S2/(S2 + kn3))
# dS3/dt
dS3 = k6 * (S4/(kn7 + S4)) \
- k3 * (S2*S3/(kn4 + S3))
# dS4/dt
dS4 = k3 * (S2*S3/(kn4 + S3)) \
- k6 * (S4/(kn7 + S4)) \
+ k5 * (S5/(kn6 + S5)) \
- k4 * (S2*S4/(kn5 + S4))
# dS5/dt
dS5 = k4 * (S2*S4/(kn5 + S4)) \
- k5 * (S5/(kn6 + S5))
# dS6/dt
dS6 = k10 * (S7/(kn11 + S7)) \
- k7 * (S5*S6/(kn8 + S6))
# dS7/dt
dS7 = k7 * (S5*S6/(kn8 + S6)) \
- k10 * (S7/(kn11 + S7)) \
- k8 * (S5*S7/(kn9 + S7)) \
+ k9 * (S8/(kn10 + S8))
# dS8/dt
dS8 = k8 * (S5*S7/(kn9 + S7)) \
- k9 * (S8/(kn10 + S8))
return np.array([dS1, dS2, dS3, dS4, dS5, dS6, dS7, dS8], dtype=float)
def solve_model_mapk_simple(initial_conditions, time_points, klin, kn, n1):
def ode(t, x):
return model_mapk_simple(x, klin, kn, n1)
sol = solve_ivp(ode,
[time_points[0], time_points[-1]],
initial_conditions,
t_eval=time_points,
method="RK45",
rtol=1e-6)
return sol
In [70]:
# klin = [k1, k2, k3, k4, k5, k6, k7, k8, k9, k10]
# kn = [kn1, kn2, kn3, kn4, kn5, kn6, kn7, kn8, kn9, kn10, kn11]
# n = n1 # only the first reaction uses a Hill exponent
# s_star_mapk = [(1,1,1,1,1,1,1,1)]
klin = [
2.5, # k1
0.25, # k2
0.025, # k3
0.025, # k4
0.75, # k5
0.75, # k6
0.025, # k7
0.025, # k8
0.5, # k9
0.5 # k10
]
kn = [
9, # kn1 = KI
10, # kn2 = K1
8, # kn3 = K2
15, # kn4
15, # kn5
15, # kn6
15, # kn7
15, # kn8
15, # kn9
15, # kn10
15 # kn11
]
n1 = 1
initial_conditions = [
100, # S1
0, # S2
300, # S3
0, # S4
0, # S5
300, # S6
0, # S7
0 # S8
]
In [71]:
sol = solve_model_mapk_simple(initial_conditions,
time_points=np.linspace(0, 5000, 5001),
klin=klin,
kn=kn,
n1=n1)
In [72]:
# Plot
plt.figure(figsize=(10, 6))
plt.plot(sol.t, sol.y.T)
plt.xlim(0, 5000)
plt.xlabel('Time')
plt.ylabel('Concentration')
plt.title('Concentration vs Time for Simplified MAPK Model')
plt.legend([f'S{i+1}' for i in range(8)], bbox_to_anchor=(1.05, 1), loc='upper left')
plt.tight_layout()
plt.show()
In [73]:
# Jacobian matrix
j_matrix_mapk_simple = numdifftools.Jacobian(
lambda y: model_mapk_simple(y, klin, kn, n1),
method="central"
)(initial_conditions)
signed_jacobian_mapk_simple = np.sign(j_matrix_mapk_simple)
print("Signed Jacobian matrix for simplified MAPK model at initial conditions:")
print(signed_jacobian_mapk_simple)
Signed Jacobian matrix for simplified MAPK model at initial conditions: [[-1. 1. 0. 0. 0. 0. 0. 1.] [ 1. -1. 0. 0. 0. 0. 0. -1.] [ 0. -1. 0. 1. 0. 0. 0. 0.] [ 0. 1. 0. -1. 1. 0. 0. 0.] [ 0. 0. 0. 0. -1. 0. 0. 0.] [ 0. 0. 0. 0. -1. 0. 1. 0.] [ 0. 0. 0. 0. 1. 0. -1. 1.] [ 0. 0. 0. 0. 0. 0. 0. -1.]]
In [75]:
# Loop detection
s_star_mapk = [initial_conditions]
res_tab_mapk_simple = ld.find_loops_vset(model_mapk_simple,
vset=s_star_mapk,
klin=klin,
kn=kn,
n1=n1,
max_num_loops=50)
print(res_tab_mapk_simple['loop_rep'][0])
loop length sign 0 (0, 1, 0) 2 1 1 (0, 0) 1 -1 2 (1, 1) 1 -1 3 (3, 3) 1 -1 4 (4, 4) 1 -1 5 (6, 6) 1 -1 6 (7, 7) 1 -1
In [85]:
type(res_tab_mapk_simple['loop_rep'][0])
Out[85]:
pandas.core.frame.DataFrame
In [86]:
# Detect durations of the fluctuations
# def compute_fluctuation_durations(loop_rep, j_matrix):
# # No library function yet, so we implement a simple version here
# return durations
# durations = compute_fluctuation_durations()
Example 6: Calcium oscillations model¶
The reactions (flows) constituting the model are given by $$ \begin{aligned} \nu_1 &= k_1,\\[6pt] \nu_2 &= k_2,\\[6pt] \nu_3 &= k_3 \cdot \frac{S_1^{n_1}}{kn_1^{\,n_1} + S_1^{n_1}},\\[10pt] \nu_4 &= k_4 \cdot \frac{S_2^{n_2}}{kn_2^{\,n_2} + S_2^{n_2}} \cdot \frac{S_1^{n_3}}{kn_3^{\,n_3} + S_1^{n_3}},\\[10pt] \nu_5 &= k_5 \cdot S_2,\\[6pt] \nu_6 &= k_6 \cdot S_1. \end{aligned} $$
The model is composed from these six reactions by $$ \begin{aligned} \frac{dS_1}{dt} &= \nu_1 + \nu_2 - \nu_3 + \nu_4 + \nu_5 - \nu_6,\\[6pt] \frac{dS_2}{dt} &= \nu_3 - \nu_4 - \nu_5. \end{aligned} $$
In [5]:
def model_calcium_oscillation(y, klin, kn, n):
S1, S2 = y
k1, k2, k3, k4, k5, k6 = klin
kn1, kn2, kn3 = kn
n1, n2, n3 = n
v1 = k1
v2 = k2
v3 = k3 * (S1**n1) / (kn1**n1 + S1**n1)
v4 = k4 * (S2**n2)/(kn2**n2 + S2**n2) * (S1**n3)/(kn3**n3 + S1**n3)
v5 = k5 * S2
v6 = k6 * S1
dS1 = v1 + v2 - v3 + v4 + v5 - v6
dS2 = v3 - v4 - v5
return np.array([dS1, dS2])
def solve_calcium_oscillation(initial_conditions, time_points, klin, kn, n):
def model(t, x):
return model_calcium_oscillation(x, klin, kn, n)
sol = solve_ivp(model, [time_points[0], time_points[-1]], initial_conditions,
t_eval=time_points, method='RK45', rtol=1e-6)
return sol
In [6]:
# Parameters and initial conditions for calcium oscillation model taken from supplement material
initial_conditions_calcium = [[0.3920, 1.6456]]
klin_calcium = [1, 7.3, 65, 500, 1, 10]
kn_calcium = [1, 2, 0.9]
n_calcium = [2, 2, 4]
time_points = np.linspace(0, 3, 31)
In [ ]:
# Parameters and initial conditions for calcium oscillation model that I varies
initial_conditions_calcium = [[0.3920, 1.6456]]
klin_calcium = [1, 7.3, 65, 500, 1, 10]
kn_calcium = [2, 2, 2]
n_calcium = [1, 2, 4]
time_points = np.linspace(0, 3, 31)
In [7]:
sol = solve_calcium_oscillation(initial_conditions=initial_conditions_calcium[0],
time_points=time_points,
klin=klin_calcium,
kn=kn_calcium,
n=n_calcium)
In [8]:
plt.figure(figsize=(8, 6))
for i, c in enumerate(colors[:2]):
plt.plot(time_points, sol.y[i], lw=2.2, color=c, label=f"S{i+1}")
plt.xlabel('Time (a.u.)', weight='bold')
plt.xlim(0, 3)
plt.xticks(np.arange(0, 3.5, 0.5))
# Show tick marks in both left and right y-axes
plt.gca().yaxis.set_ticks_position('both')
plt.ylabel('Concentration (a.u.)', weight='bold')
# plt.yscale('log')
# plt.title(f'Parameters: {params}')
# Make legend outside the plot
plt.legend(loc='upper left', bbox_to_anchor=(1, 1))
# plt.suptitle('Figure 2D Variations: Concentration vs Time for Old Model with Loops (fixed)')
plt.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.title('Calcium Oscillation Model: Concentration vs Time (original params)', weight='bold')
#Show in 300 dpi
plt.gcf().set_dpi(300)
plt.show()
In [13]:
# initial_conditions_calcium = [[0.3920*10**(-6), 1.6456*10**(-6)]]
# klin_calcium = [1*10**(-6), 7.3*10**(-6), 65*10**(-6), 500*10**(-6), 1, 10]
# kn_calcium = [1*10**(-6), 2*10**(-6), 0.9*10**(-6)]
# n_calcium = [2, 2, 4]
# time_points = np.linspace(0, 3, 31)
initial_conditions_calcium = [[0.3920, 1.6456]]
klin_calcium = [1, 7.3, 65, 500, 1, 10]
kn_calcium = [2, 2, 2]
n_calcium = [1, 2, 4]
time_points = np.linspace(0, 10, 101)
In [14]:
sol = solve_calcium_oscillation(initial_conditions=initial_conditions_calcium[0],
time_points=time_points,
klin=klin_calcium,
kn=kn_calcium,
n=n_calcium)
In [15]:
plt.figure(figsize=(8, 6))
for i, c in enumerate(colors[:2]):
plt.plot(time_points, sol.y[i], lw=2.2, color=c, label=f"S{i+1}")
plt.xlabel('Time (a.u.)', weight='bold')
plt.xlim(0, 10)
plt.xticks(np.arange(0, 11, 1))
# Show tick marks in both left and right y-axes
plt.gca().yaxis.set_ticks_position('both')
plt.ylabel('Concentration (a.u.)', weight='bold')
# plt.yscale('log')
# plt.title(f'Parameters: {params}')
# Make legend outside the plot
plt.legend(loc='upper left', bbox_to_anchor=(1, 1))
plt.title('Calcium Oscillation Model: Concentration vs Time (modified params)', weight='bold')
plt.tight_layout(rect=[0, 0.03, 1, 0.95])
#Show in 300 dpi
plt.gcf().set_dpi(300)
plt.show()
In [21]:
# vset=[(1.0, 1.0)],
# klin=(0.5,0.3,1.0,1.5,0.2,0.1),
# kn=(0.5,0.5,0.5),
# n=(2,2,2)
res_tab = ld.find_loops_vset(model_calcium_oscillation,
vset=initial_conditions_calcium,
klin=klin_calcium,
kn=kn_calcium,
n=n_calcium,
max_num_loops=50)
# The loop list, a pandas dataframe, is accessed like this:
res_tab['loop_rep'][0]
Out[21]:
| loop | length | sign | |
|---|---|---|---|
| 0 | (0, 1, 0) | 2 | 1 |
| 1 | (0, 0) | 1 | -1 |
| 2 | (1, 1) | 1 | -1 |
In [11]:
# Jacobian matrix at the point of interest
j_matrix_calcium = numdifftools.Jacobian(model_calcium_oscillation,method="central")(
initial_conditions_calcium[0],
klin=klin_calcium,
kn=kn_calcium,
n=n_calcium
).real
signed_jacobian_calcium = np.sign(j_matrix_calcium)
print("Signed Jacobian matrix for calcium oscillation model at point of interest:")
print(signed_jacobian_calcium)
Signed Jacobian matrix for calcium oscillation model at point of interest: [[ 1. 1.] [-1. -1.]]
In [36]:
import numpy as np
from scipy.optimize import root
def f_rhs(x, klin, kn, n):
return model_calcium_oscillation(x, klin, kn, n)
def find_equilibrium(x0, klin, kn, n):
sol = root(lambda x: f_rhs(x, klin, kn, n), x0, method="hybr")
if not sol.success:
raise RuntimeError(f"Equilibrium root-find failed: {sol.message}")
return sol.x
In [37]:
def numerical_jacobian(fun, x, eps=1e-7):
x = np.asarray(x, dtype=float)
f0 = np.asarray(fun(x))
J = np.zeros((len(x), len(x)))
for j in range(len(x)):
xp = x.copy(); xp[j] += eps
xm = x.copy(); xm[j] -= eps
fp = np.asarray(fun(xp))
fm = np.asarray(fun(xm))
J[:, j] = (fp - fm) / (2*eps)
return J
def equilibrium_stability(x_star, klin, kn, n):
J = numerical_jacobian(lambda x: f_rhs(x, klin, kn, n), x_star)
eigvals = np.linalg.eigvals(J)
return J, eigvals
In [45]:
x_star_calcium = find_equilibrium(initial_conditions_calcium[0],
klin=klin_calcium,
kn=kn_calcium,
n=n_calcium)
print("Equilibrium point (calcium model):", x_star_calcium)
Equilibrium point (calcium model): [0.83 0.74820287]
In [46]:
j_matrix_calcium_eq = numerical_jacobian(lambda x: f_rhs(x, klin_calcium, kn_calcium, n_calcium), x_star_calcium)
signed_jacobian_calcium_eq = np.sign(j_matrix_calcium_eq)
print("Jacobian matrix at equilibrium point (calcium model):")
print(j_matrix_calcium_eq)
print("Signed Jacobian matrix for calcium oscillation model at equilibrium point:")
print(signed_jacobian_calcium_eq)
Jacobian matrix at equilibrium point (calcium model): [[ 24.22375257 61.41680616] [-34.22375256 -61.41680615]] Signed Jacobian matrix for calcium oscillation model at equilibrium point: [[ 1. 1.] [-1. -1.]]
In [47]:
j_matrix_calcium_eq, eigvals_calcium_eq = equilibrium_stability(x_star_calcium, klin=klin_calcium, kn=kn_calcium, n=n_calcium)
print("Eigenvalues at equilibrium point (calcium model):", eigvals_calcium_eq)
Eigenvalues at equilibrium point (calcium model): [-18.59652679+16.38100279j -18.59652679-16.38100279j]
In [48]:
from scipy.integrate import solve_ivp
from scipy.signal import find_peaks
def solve_calcium_oscillation_dense(x0, t_end, dt, klin, kn, n):
t_eval = np.arange(0, t_end + dt, dt)
def model(t, x):
return model_calcium_oscillation(x, klin, kn, n)
sol = solve_ivp(model, [t_eval[0], t_eval[-1]], x0,
t_eval=t_eval, method='RK45', rtol=1e-7, atol=1e-9)
return sol
def detect_limit_cycle_from_timeseries(t, x, transient_frac=0.5, var_index=0,
min_peaks=6, amp_tol=1e-4):
"""
Returns: (is_cycle, amplitude, period_est)
"""
N = len(t)
cut = int(transient_frac * N)
tt = t[cut:]
xx = x[var_index, cut:]
# Peak detection (you can tune `distance` if needed)
peaks, _ = find_peaks(xx)
if len(peaks) < min_peaks:
return False, 0.0, np.nan
# Estimate period from last few peaks
peak_times = tt[peaks]
dT = np.diff(peak_times)
if len(dT) < (min_peaks - 1):
return False, 0.0, np.nan
period_est = np.median(dT[-(min_peaks-2):])
# Amplitude estimate (peak-to-trough)
amp = float(np.max(xx) - np.min(xx))
if amp < amp_tol:
return False, amp, period_est
# Extra check: variability persists (not decaying)
# Compare early vs late in the "post-transient" window
mid = len(xx)//2
amp1 = np.max(xx[:mid]) - np.min(xx[:mid])
amp2 = np.max(xx[mid:]) - np.min(xx[mid:])
if amp2 < 0.5 * amp1: # suggests decay
return False, amp, period_est
return True, amp, period_est
In [49]:
def sweep_one_parameter(x0, klin, kn, n, which=("klin", 1),
values=np.linspace(0.5, 20, 80),
t_end=200, dt=0.02, var_index=0):
results = []
for val in values:
klin2 = list(klin); kn2 = list(kn); n2 = list(n)
group, idx = which
if group == "klin":
klin2[idx] = float(val)
elif group == "kn":
kn2[idx] = float(val)
elif group == "n":
n2[idx] = float(val)
else:
raise ValueError("which must be ('klin',i), ('kn',i), or ('n',i)")
sol = solve_calcium_oscillation_dense(x0, t_end, dt, klin2, kn2, n2)
t = sol.t
X = sol.y
is_cycle, amp, period = detect_limit_cycle_from_timeseries(t, X, var_index=var_index)
# extrema after transient for diagram
cut = int(0.5 * len(t))
post = X[var_index, cut:]
x_min = float(np.min(post))
x_max = float(np.max(post))
results.append({
"param": float(val),
"is_cycle": bool(is_cycle),
"amp": float(amp),
"period": float(period),
"x_min": x_min,
"x_max": x_max,
})
return results
In [50]:
sweep_one_parameter_results = sweep_one_parameter(
x0=initial_conditions_calcium[0],
klin=klin_calcium,
kn=kn_calcium,
n=n_calcium,
which=("kn", 2),
values=np.linspace(0.5, 20, 80),
t_end=200,
dt=0.02,
var_index=0
)
In [51]:
sweep_one_parameter_results[:5] # Show first 5 results
Out[51]:
[{'param': 0.5,
'is_cycle': False,
'amp': 1.2874821253561208e-07,
'period': 0.0799999999999983,
'x_min': 0.8299999270052814,
'x_max': 0.8300000557534939},
{'param': 0.7468354430379747,
'is_cycle': False,
'amp': 1.6121071944041887e-07,
'period': 0.04999999999999716,
'x_min': 0.8299999055148357,
'x_max': 0.8300000667255552},
{'param': 0.9936708860759493,
'is_cycle': False,
'amp': 1.8390659028977296e-07,
'period': 0.21999999999999886,
'x_min': 0.8299999200077249,
'x_max': 0.8300001039143152},
{'param': 1.240506329113924,
'is_cycle': True,
'amp': 1.870155214255138,
'period': 0.6200000000000045,
'x_min': 0.3996457383673823,
'x_max': 2.2698009526225205},
{'param': 1.4873417721518987,
'is_cycle': True,
'amp': 3.8780667000241578,
'period': 1.25,
'x_min': 0.3723368455513571,
'x_max': 4.250403545575515}]
In [52]:
# --- Convert list-of-dicts to arrays
params = np.array([r["param"] for r in sweep_one_parameter_results], dtype=float)
x_min = np.array([r["x_min"] for r in sweep_one_parameter_results], dtype=float)
x_max = np.array([r["x_max"] for r in sweep_one_parameter_results], dtype=float)
amp = np.array([r["amp"] for r in sweep_one_parameter_results], dtype=float)
period = np.array([r["period"] for r in sweep_one_parameter_results], dtype=float)
is_cyc = np.array([r["is_cycle"] for r in sweep_one_parameter_results], dtype=bool)
# --- 1) Bifurcation-style plot: extrema vs parameter
plt.figure()
plt.plot(params, x_min, marker=".", linestyle="None", label="post-transient min")
plt.plot(params, x_max, marker=".", linestyle="None", label="post-transient max")
# Optional: visually distinguish points where a cycle was detected
plt.plot(params[~is_cyc], x_min[~is_cyc], marker="x", linestyle="None", label="no-cycle (min)")
plt.plot(params[~is_cyc], x_max[~is_cyc], marker="x", linestyle="None", label="no-cycle (max)")
plt.xlabel("Parameter value (kn[2])")
plt.ylabel("State variable S1 (var_index=0) extrema")
plt.title("Bifurcation-style diagram (extrema after transient)")
plt.legend()
plt.tight_layout()
plt.show()
# --- 2) Amplitude vs parameter
plt.figure()
plt.plot(params, amp, marker="o", linestyle="-")
plt.xlabel("Parameter value (kn[2])")
plt.ylabel("Amplitude (max-min, post-transient)")
plt.title("Oscillation amplitude vs parameter")
plt.tight_layout()
plt.show()
# --- 3) Period vs parameter (only where a cycle is detected)
plt.figure()
plt.plot(params[is_cyc], period[is_cyc], marker="o", linestyle="None")
plt.xlabel("Parameter value (kn[2])")
plt.ylabel("Estimated period")
plt.title("Limit cycle period vs parameter (cycle points only)")
plt.tight_layout()
plt.show()
In [55]:
sweep_one_parameter_results = sweep_one_parameter(
x0=initial_conditions_calcium[0],
klin=klin_calcium,
kn=kn_calcium,
n=n_calcium,
which=("kn", 2),
values=np.linspace(0, 3, 61),
t_end=200,
dt=0.02,
var_index=0
)
In [54]:
np.linspace(0, 3, 61)
Out[54]:
array([0. , 0.05, 0.1 , 0.15, 0.2 , 0.25, 0.3 , 0.35, 0.4 , 0.45, 0.5 ,
0.55, 0.6 , 0.65, 0.7 , 0.75, 0.8 , 0.85, 0.9 , 0.95, 1. , 1.05,
1.1 , 1.15, 1.2 , 1.25, 1.3 , 1.35, 1.4 , 1.45, 1.5 , 1.55, 1.6 ,
1.65, 1.7 , 1.75, 1.8 , 1.85, 1.9 , 1.95, 2. , 2.05, 2.1 , 2.15,
2.2 , 2.25, 2.3 , 2.35, 2.4 , 2.45, 2.5 , 2.55, 2.6 , 2.65, 2.7 ,
2.75, 2.8 , 2.85, 2.9 , 2.95, 3. ])
In [56]:
# --- Convert list-of-dicts to arrays
params = np.array([r["param"] for r in sweep_one_parameter_results], dtype=float)
x_min = np.array([r["x_min"] for r in sweep_one_parameter_results], dtype=float)
x_max = np.array([r["x_max"] for r in sweep_one_parameter_results], dtype=float)
amp = np.array([r["amp"] for r in sweep_one_parameter_results], dtype=float)
period = np.array([r["period"] for r in sweep_one_parameter_results], dtype=float)
is_cyc = np.array([r["is_cycle"] for r in sweep_one_parameter_results], dtype=bool)
# --- 1) Bifurcation-style plot: extrema vs parameter
plt.figure()
plt.plot(params, x_min, marker=".", linestyle="None", label="post-transient min")
plt.plot(params, x_max, marker=".", linestyle="None", label="post-transient max")
# Optional: visually distinguish points where a cycle was detected
plt.plot(params[~is_cyc], x_min[~is_cyc], marker="x", linestyle="None", label="no-cycle (min)")
plt.plot(params[~is_cyc], x_max[~is_cyc], marker="x", linestyle="None", label="no-cycle (max)")
plt.xlabel("Parameter value (kn[2])")
plt.ylabel("State variable S1 (var_index=0) extrema")
plt.title("Bifurcation-style diagram (extrema after transient)")
plt.legend()
plt.tight_layout()
plt.show()
# --- 2) Amplitude vs parameter
plt.figure()
plt.plot(params, amp, marker="o", linestyle="-")
plt.xlabel("Parameter value (kn[2])")
plt.ylabel("Amplitude (max-min, post-transient)")
plt.title("Oscillation amplitude vs parameter")
plt.tight_layout()
plt.show()
# --- 3) Period vs parameter (only where a cycle is detected)
plt.figure()
plt.plot(params[is_cyc], period[is_cyc], marker="o", linestyle="None")
plt.xlabel("Parameter value (kn[2])")
plt.ylabel("Estimated period")
plt.title("Limit cycle period vs parameter (cycle points only)")
plt.tight_layout()
plt.show()
Additional¶
Exporting environment details to yaml file¶
In [4]:
! conda env export > environment.yml
In [17]:
# --- Export notebook to HTML (recommended) ----------------------------------
import subprocess, sys, pathlib
NB_NAME = "test.ipynb" # <— set once to your actual filename
OUT_DIR = pathlib.Path("./")
# OUT_DIR = OUT_DIR.resolve()
# OUT_DIR.mkdir(parents=True, exist_ok=True)
subprocess.run([
sys.executable, "-m", "jupyter", "nbconvert",
"--to", "html",
"--output-dir", str(OUT_DIR),
NB_NAME
], check=False)
print("HTML saved to", OUT_DIR / (pathlib.Path(NB_NAME).stem + ".html"))
HTML saved to test.html